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We present a detailed study of the effect of time delay on the collective dynamics of coupled limit 
cycle oscillators at Hopf bifurcation. For a simple model consisting of just two oscillators with a 
time delayed coupling, the bifurcation diagram obtained by numerical and analytical solutions shows 
significant changes in the stability boundaries of the amplitude death, phase locked and incoherent 
regions. A novel result is the occurrence of amplitude death even in the absence of a frequency 
mismatch between the two oscillators. Similar results are obtained for an array of N oscillators with 
a delayed mean field coupling and the regions of such amplitude death in the parameter space of 
the coupling strength and time delay are quantified. Some general analytic results for the N — > oo 
(thermodynamic) limit are also obtained and the implications of the time delay effects for physical 
applications are discussed. 
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I. INTRODUCTION 



The collective behavior of a large assembly of coupled nonlinear oscillators provides valuable clues for understanding 
the complex dynamics of systems with many degrees of freedom. This has been one of the major motivations for the 
recent large scale interest both experimentally [Q and numerically j2j in the study of such simple mathematical models 
and their application to a wide variety of physical and biological problems. Examples of such applications include 
interactions of arrays of Josephson junctions GLB], semiconductor lasers fsjH, charge density waves (7|, phase locking of 
relativistic magnetrons Bclousov-Zhabotinskii reactions in coupled Brusselator models P~p^| and neural oscillator 
models for circadian pacemakers fl3|| . One of the prominent cooperative phenomena, that was first highlighted by 
Winfrcc |l4| in a simple model of weakly coupled limit cycle oscillators, is that of frequency entrainment or synchro- 
nization of the diverse frequencies of the oscillator assembly to a single common frequency Jl5|,|l6| . This happens in 
a spontaneous and abrupt fashion once the coupling strength exceeds a critical threshold. Real life examples of such 
behavior abound in nature e.g. the synchronous flashing of a swarm of fire flies, the chirping of crickets in unison or 
the electrical synchrony in cardiac cells. When the coupling between the oscillators is comparable to the attraction 
to their individual limit cycles, other interesting phenomena can occur JT7[-pO| which involve the amplitudes of the 
individual oscillators. For example, if the coupling is sufficiently strong and the spread in the natural frequencies of 
the oscillators sufficiently broad, the oscillators can suffer an amplitude quenching or death ]2l|-p3||. Such behavior 
has been observed in experiments of coupled chemical oscillator systems e.g. coupled Belousov-Zhabotinskii reactions 
carried out in coupled tank reactors p4[ ]. Other collective phenomena that these coupled oscillator models display 
include partial synchronization, phase trapping, large amplitude Hopf oscillations and even chaotic behavior |25| , p6[| - 
all of which have been discussed widely in the literature. 

The goal of our present study is to examine the effect of time delay on the collective dynamics of coupled oscil- 
lator systems. Time delay is ubiquitous in most physical and biological systems like optical bistable devices [p7[, 



electromechanical systems |28 , predator-prey models pty , and physiological systems |5C|[3l|. They can arise from 
finite propagation speeds of signals, for example, or from finite processing times in synapses, finite reaction times in 
chemical processes and so on. Surprisingly most past studies on coupled oscillator systems have not considered the 
effect of time delay. The work of Schuster and Wagner |32| , Niebur et al. |?3| , Nakamura, et al. |34|] and Kim, et al. 
p5| , are the only ones, that we are aware of, where they have tried to incorporate time delay effects in the context of 
the coupled oscillator problem. However, they have restricted themselves to the simplest of models, that of coupled 
limit cycle oscillators in the weak coupling limit where only the phase information is retained and phenomena like 
amplitude death cannot occur. Our aim is to extend this study into the strong coupling regime where both phase 
and amplitude responses need to be retained and to investigate the effect of time delay on the various collective 
responses of such a model system. To keep the analysis simple we begin with just two limit cycle oscillators that 
are coupled with a time delay. For such a model the bifurcation diagram is easy to obtain both analytically as well 
as numerically. Furthermore it allows us a detailed comparison with the past work of Aronson et al. who have 
analyzed an analogous model but without any time delay in the coupling. After identifying and briefly discussing the 



1 



results of time delay effects in this simple N — 2 model we next proceed to analyze a large assembly of coupled limit 
cycle oscillators. For this study, we construct a model which is a generalization of the mean field model of Matthews 
et al. [ p5|j2^ ] and where the coupling term is suitably modified to introduce time delay. This model is explored in 
detail by extensive numerical solutions and linear stability analysis around fixed points. The limit of N — > oo is 
particularly interesting and permits some explicit analytic results. Our principal focus is on cooperative phenomena 
like amplitude death and frequency locking and we find that both these states are significantly influenced by time 
delay effects. One of the surprising and dramatic results is that in the presence of time delay amplitude death can 
occur even for zero frequency mismatch between the oscillators (i.e. for identical oscillators). This is in sharp contrast 
to the situation with no time delay where all previous numerical and analytical works have emphasized the need to 
have a broad frequency spread for amplitude death to occur. A brief report of this result has been published by us 
elsewhere |3t| . In this paper we give a more detailed and complete description of this phenomenon. We also report 
on other newer findings related to time delay induced effects in the collective regimes corresponding to phase locked 
and chaotic states. 



The organization of the paper is as follows. In the next section (Section ||) we analyze the model of two limit 
cycle oscillators that have a time delayed coupling and compare and contrast our results with the previous work of 
Aronson et al. In Section [II , we describe the more general N coupled oscillator model and present numerical as well 
as analytic results for the different collective states. This includes the amplitude death region, the phase locked region 
and the so called chaotic regime. Some explicit results for the N — ► oo (thermodynamic) limit are also presented. 
Section IV provides a summary of our results and a brief discussion on possible future extensions of this work. 



II. TWO DELAY COUPLED OSCILLATORS 



A. Model Equations 



For the fundamental oscillator unit of our model we choose the simple limit cycle oscillator described by the equation, 

Z j (t) = (l + iu j -\Z i {t)\ 3 )Z j (t) (I) 

where Zj is the complex amplitude of the jth oscillator. Each oscillator has a stable limit cycle of unit amplitude 
| Zj | =1 with angular frequency u)j . We consider a simple model in which two of them arc coupled linearly to each 
other as follows, 

Zi(t) = (1 + twi- | Zt(t) | 2 ) Z 1 (t) + K[Z 2 {t~r) - Z x {t)\, (2a) 
Z 2 (t) = (l + iw 2 - | Z 2 (t) | 2 ) Z 2 (t)+ K [Zi(t-r) - Z 2 (t)], (2b) 

where K > is the coupling strength and r > is a measure of the time delay. These model equations are a 
direct generalization of the set of equations studied by Aronson, et al. who do not have any time delay. The 
coupled equations represent the interaction between weakly nonlinear oscillators (that are near a supercritical Hopf 
bifurcation) and whose coupling strength is comparable to the attraction of the limit cycles. It _Js important then 
to retain both the phase and amplitude response of the oscillators. Writing Zj = rje l6j , Eqs. 
expressed in polar form as, 



h = n{l - K - r{) + Kr 2 {t - t) cos[6 2 {t - t) - Q x ], (3a) 

r 2 = r 2 (l - K - r\) + Kn(t - r) cos[9x{t - r) - 2 ], (3b) 

fli = h>i + K T2{t - T) sm[d 2 (t - r) - (3c) 

n 

9 2 = lo 2 + K ri ^~ T > sin[0i (t - t) - 6 2 ]. (3d) 
r 2 



For t = 0, the work of Aronson et al. [|lj shows that the nonlinear Eqs.( pa] - |3d|) have a variety of stationary and 
non-stationary solutions which depend on the strength of the coupling parameter K and the frequency mismatch 
between the oscillators A =| wi — u> 2 |. For extremely weak coupling (K — * 0) and large A, the oscillators behave 
independently and the long term behavior is a nonstationary incoherent state in which the relative phase of the two 
oscillators moves through all phases. Such a state is also called a phase drift state. With increasing coupling strength 
two important classes of stationary solutions are possible. One of them is amplitude death in which the oscillators 



.2b) can also be 
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pull each other off their limit cycles and collapse into the origin (n = r 2 = 0) as t — > oo. The other collective state 
is called frequency locking or mutual entrainment in which the two oscillators synchronize to a common frequency 
and the time asymptotic state is one of coherent or collective oscillation. The distribution of these solutions can be 
neatly represented in a phase diagram (bifurcation diagram) in the A — K space. Fig. 1, reproduced from the work 
of Aronson et al. summarizes the above discussion. Region I represents the amplitude death region, region III is 
the phase drift regime and region II marks the phase locked state. We now analyze Eqs.(^a -3d) for finite values of r 
and examine the effect of r on the conditions for the onset of these states and changes if any in the basic properties 
of these states. In the following subsections we discuss the results for the amplitude death and phase locked solutions. 



B. Amplitude Death 



As is clear from Eqs.( pa| , pb[ ), the origin Zj = (J = 1,2) is always a fixed point of the system. The question to 
consider is whether this is a stable fixed point in which case all amplitudes of the oscillators would collapse into the 
origin as t — > oo. In the absence of time delay this state occurs in the region K > 1 for A > 2. To determine the 
onset conditions of this state in the presence of time delay (r ^ 0) we linearize Eqs.(^a|,|l]) around Zj = to obtain 
the characteristic equation, 



det{A — XI) = 



(4) 



where A, the linearized matrix of the Eqs.(|2a|) and (2b), is given by 



A = 



iu>i Ke 



-At 



Ke 



-At 



a + iuj2 



(5) 



I is the identity matrix, a = 1 — K and the perturbations are assumed to have a time dependence proportional to e At . 
Eqn.(0) can be written in the form 



(a — A + iu>i)(a — A + itoz) — K 2 e 



-2Ar 



0. 



(6a) 



or 



A 2 - 2(a + iCo)\ + (h + ib 2 ) + ce~ 2XT = (6b) 

where b\ = a 2 — Co 2 + A 2 /4, b 2 = 2au>, A =| lo% — uo 2 \, ui — (uii + u>2)/2 and c = —K 2 . This is a transcendental 
equation having an infinite number of roots and we wish to study the movement of the eigenvalues in the parametric 
plane of (K, A) and (K, r). Setting A = a + if3, where a and (3 are real, the amplitude death region corresponds to 
the region in which a < 0. The marginal stability curves or the critical curves are thus obtained by requiring that 
a = 0, i.e. A = if]. Substituting in (|6a|), the equations defining the critical curves are thus given by, 

(/? - uj) 2 - A 2 /4 -a 2 +K 2 cos(2/3r) = 0, (7a) 
2a{/3 -u))~K 2 sin(2/3r) = 0. (7b) 

We first briefly describe the case of r = 0, in order to appreciate the changes brought about by finite time de- 
lay. Setting r = in (^) we obtain the conditions K = 1 and (3 = uj. Substituting for (3 in ( [7a] ) we obtain 

K = 7(A) = + A-). So the critical curves in this case are K = 1 and K = 7(A) in agreement with the work of 
Aronson, et al. pl| and as illustrated in Fig. 1. It is appropriate to distinguish between two regions in the (if, A) 
space, namely (i) A > 2 and (ii) A < 2. When A > 2 the stable region of the origin (amplitude death region) 
is bounded by K = 1 and K = 7(A). The eigenvalues in this particular case can be written down (from ([3l])) as 
A = 1 — K ± yj K 2 — A 2 /4 ± iu). On the boundary K = 1, the origin loses stability in a Hopf bifurcation. Two pairs 
of eigenvalues cross into the right half plane. On the boundary K = 7(A) a pair of eigenvalues crosses into the right 
hand side of the complex eigenvalue plane, giving rise to a single frequency, which corresponds to the phase locked 
state of the system. In the second case, when A < 2 there is no amplitude death. However the critical curves give 
a boundary on which an unstable fixed point is born as one moves to the left of K = 7(A) curve. Note also that 
the boundaries of the three regions meet in a highly degenerate manner in the single point K — 1, A = 2. Another 
distinguishing feature of the r = case is that the critical curves are independent of the mean frequency Co. In fact 
one can set Co — (which is equivalent to transforming to a frame rotating at the mean frequency) and carry out 
the same analysis without any loss of generality. This property follows from the original symmetry of the coupled 
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equations. 



When t ^ 0, this symmetry is lost and the critical curves are no longer independent of Q as seen from Eqs.(|7a|, ffb| ). 
In comparing the phase diagrams of the r ^ case with that of the Aronson et al. [M diagram we therefore need to 
always mention the specific value of the mean frequency parameter. We now briefly describe our numerical procedure 
for solving Eqs. (|7a|, 7b) in order to plot graphically the critical curves in the {K, A) space. To eliminate (3 between 
the two equations it is more convenient to write the equations in the following parametric form, 

F = ([3 — Q)/ sin(2/3r), (8a) 



K = K± = —F ± \Jf 2 + 2F, (8b) 
A 2 = -4a 2 + 4(/3 - uf + AK 2 cos(2/3r). (8c) 



where (|8b|) corresponds to the two roots of K for the quadratic equation (7b). Equations ( pa| - |8c| ) represent two sets 
of curves arising due to the ± signs. Let us denote the set of curves arising due to the "+" sign by S+ = S+(K + , A) 
and the curves due to the "— " sign by 5- = S-(K~, A). The curves S+ and S- are obtained by choosing an 
interval of (3 and correspondingly evaluating K and A and thus eliminating (3. The function F has singularities at 
(3 = (3 n = mr/ (2r) where n is an integer. Each interval / = (f3 n , (3 n +i) provides a part of the phase curves in (K, A) 
plane. In the (K, A) space, the curves S- exist between K = 1 and K = 2 and the curves S+ exist outside of this 
region in the intervals 1/2 < K < 1 and K > 2. And at higher values of r, 5_ produces curves which could intersect. 
For small values of the parameter r close to 0, the amplitude death region is bounded by the curve S- when K < 1 
and S+ when K > 1. At higher values of r the boundary of the death region falls below K = 1 in which region the 
boundary is specified by the curves S+. The transcendental equations ( pafpq ) must be studied for each parametric 
value of r since the curves in (K, A) space become more complicated as the parameter r is increased. 

We now present our results of the amplitude death region for Co — 10 for various values of r as obtained by the 
numerical prescription described above. From the series of diagrams in Fig. 2 we notice that with the introduction 
of finite r the point (K = 1, A = 2) no longer has a degenerate character and the critical curves begin separating and 
distorting in a continuous manner. The amplitude death region grows in size as the value of r is increased from and 
the curve S+ defined by the interval of (3 € (0, Co) starts bending down below the A = 2\J2K — 1 curve towards the 
A = axis. For a critical value of r = t c it touches the A = axis and the region of death on the axis lies between 
the two points of intersection K\ and K2 of S+ with A = 0. This death region has a finite width K2 — K\ for a 
range of values of the delay parameter t. This phenomenon of death of identical oscillators is a novel result purely 
induced by the temporal delay in the coupling of the oscillators. This behavior persists for a range of t after which 
the bifurcation curve lifts up from the A = line and starts moving upward. 

To quantitatively study this region of amplitude death for identical oscillators, we take a more detailed look at the 
trajectory of the two bounding points K\ and K2 in the parametric space of (K, r) for A = 0. Let this trajectory be 
called Tb(K) for which we require that all the eigenvalues of the original transcendental equation lie in the left half 
plane of the complex eigenvalue plane when r > r c . To identify such a curve from all the permissible multiple critical 
curves in the parametric space (K,t), it is simpler to start with the original equation (|6a| ) and set lo\ = U2 = in 
it. The eigenvalue equation now simplifies to: 

A = 1 - K + iu ± Ke- XT (9) 

Let X — a + if3, where a and [3 are real. We assume without loss of generality that (3 > 0. In order to obtain the 
critical curves we set a — and consider both the equations arising out of " + " and " — " signs in (||) . After some 
straightforward algebra, involving the choice of the correct signs in the inversion of the cosine function, we obtain the 
following two sets of critical curves. 

I ZS\ H7T + COS -1 (1 — 1/ K) 

n = n{n,K) = , (10a) 

uj — \J2K — 1 

nsrM = iji±lh^^m. (10b) 

ui + v 2iv — 1 

where n = 0, 1, 00. We further need to know the nature of the transition of pairs of eigenvalues as it crosses these 
curves. For this it is necessary to evaluate da/dr on each of these curves and examine its sign. Setting A = a + i(3 in 
equation (m and differentiating with respect to t, it is straightforward to get, 
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(11) 



where c\ = [(1 ± Kt) 2 + (Kt sin(/3r)) 2 ] 1 , which is a real positive constant and (3 = (3± — ujq ± \J2K — 1. From the 
above equation it is easily seen that 

da r < on n if K < f(u ) 

— |o=o < > on n if K > f(uo) (12) 
" I > on T2 

where f(u>o) = (1 + u>q)/2. Thus on a ri curve, a pair of eigenvalues transits to the left half plane provided the 
coupling strength is smaller than /(luq) and to the right side if the coupling strength is greater than /(^o)- On a 
T2 branch of the critical curves however a pair of eigenvalues always crosses into the right half plane of the complex 
plane. Thus for a finite region of amplitude death to exist in the K — t plane it needs to be bounded by appropriate 
branches of the t\ and r 2 curves and condition K < f(tOo) should hold. For K > /(loq), there would be no amplitude 
death region at all. 



In Fig. 3, we illustrate the above arguments more graphically by plotting the critical curves Ti(n, K) (solid lines) 
and T2{n,K) (dashed lines) for the first few values of n and for w = 10 in the (K — r) space. It is possible to 
identify t\(0,K) as Tb(K). This curve forms the first boundary of the amplitude death region of identical oscillators 
in (K, t) space. The stability of the origin can be lost when a pair of eigenvalues makes a transition to the right half 
plane, which in the present case will occur on T2(0, K) (as can be verified from ([l2"l)). So the region enclosed by the 
intersection of Ti(0, K) and T2(0, K) forms a region of amplitude death in the K — t space which we label as a death 
island. Physically such an island represents a region in phase space where for a given value of uj and at a fixed K we 
move (by varying r) from an unstable region (corresponding to phase locked states) into a stable region as we cross 
the left boundary of the island to emerge again into a unstable region as we cross the right boundary of the island. Is 
it possible to have more than one island for a given value of ojq? To answer this question we see from Fig. 3, that the 
next curve along the r axis is T2(1,K) on which another pair of eigenvalues will make a transition to the right half 
plane. Thus when one gets to ri(l, K), which is the next in the sequence and on which a pair of eigenvalues crosses 
to the left half plane, there is still a pair left in the right half plane. Thus the region between n(l,K) and T2(2,K) 
does not constitute a death island. In general the frequency of r 2 curves is higher than n, i.e., 

ST 2 {n,K) < Sn{n,K) (13) 

where <5ri(n, K) — Ti(n + 1,K) — T\(n,K) and St2(ji,K) — T2(n + l,K) — T2(n,K). Thus the death island region 
is usually singly connected and there are no higher order islands. However the ordering of the curves depends on 
the magnitude of ujq. In the present case we do not see any multiple islands in the range 4.812 < ui < 14.438. For 
ujq > 14.438 we do see the appearance of higher order islands as shown in Fig. 4 for ujo = 30. The size of the primary 
death island is found to be a function of ujq and as we shall see in the next section it also depends on TV the number 
of coupled oscillators. Fig. 5, displays the island sizes for different values of luq. The size of the island decreases with 
decreasing frequency and vanishes below a certain threshold. 

All the above features of the amplitude death phenomenon have also been confirmed by a direct numerical solution 
of the coupled oscillator equations and excellent agreement with the analytic results have been found. 



C. Frequency Locking 

The frequency (or phase) locked solutions of the system (^a -3d) are characterized by the property that the relative 
phase of the two oscillators is a constant. The phase locked state can be described by the ansatz 9\ : 2{t) = Ot ± a/2 
where a, the phase difference between the two oscillators, and Q, the common frequency of the two oscillators 
are real constants. Substitution of this ansatz in Eqs.(|3"c|) and (|3d|) further shows that the amplitudes of the 
limit cycles remain constant in this case. Thus the phase locked solutions can be described by the representa- 
tion, ( ri ( t),r 2 (t),9i(t),92(t)) — (Ri, R 2 , fit — a/2, ft t + a/2) where are constants. Substituting such a form 



in (3a -3d), we obtain the following set of four equations from which the values of i?i,i?2>^ and a can be evaluated. 



(1 - K - RDRx + KR 2 cos(a - fir) = 0, (14a) 

(1 - K - Rl)R. 2 + KR X cos(a + VLt) = 0, (14b) 

fl = uj 1 +K— sin(a-nr), (14c) 
Ri 
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fl = u 2 - sin(a + Or). (14d) 
R2 

These equations can also be rearranged in the following form which is slightly more convenient for numerical solutions. 

Rl = 1- K + Kf!Cos(a- Or), (15a) 
R% = 1 - K + Kf 2 cos(a + QT), (15b) 

a = sin" 1 ^/sin 2 (fir) - (0 - wi)(0 - oj 2 )/K 2 , (15c) 

#2 - /i^i, (15d) 

where /j. = (fl — ui\)/(K sm{a — Or)) and /2 = — (f2 — u>2)/(K sin(a + Or)). This is a set of transcendental equa- 
tions whose solu tions descri be t he phase locked equilibria. We again first briefly discuss the r = case. Putting 
t = in Eqs. (HI) and ( |l5b|) , we see that there are two possible equilibrium solutions which are given by (i) 



R\ = R\ = p 2 which is called the symmetric equilibrium, and (ii) R\ + R\ = 1 — K, the asymmetric case. Both these 
equilibria have been studied in detail by Aronson, et al. [ pl| . The symmetric phase locked equilibria are given by 
p\ = 1-K± V / K 2 - A 2 /4, fl = (u; 1 +lu 2 )/2 and a±, where a+ = sinT 1 (A/(2K)), and a_ = vr - sin" 1 (A/(2if). Of 
these two symmetric equilibria, the set given by (p+, f2, a + ) is found to be stable, whereas the solution (p 2 ., 0, a_) is 
unstable. The asymmetric phase locked solutions turn out to be unstable. Thus for r = the only stable equilibrium 
is the one where the two oscillators are synchronized to the mean frequency and have identical amplitudes. Note that 
the amplitudes are lowered from the unity value of the uncoupled case by the amount K — J K 2 — A 2 /4. 



With fin ite t i me delay (r 7^ 0) there is a richer fare of equilibria. This is evident from the full set of transcendental 
equations (15a-15d) which permit a large number of solutions. Although it is difficult to obtain analytic forms for 
these solutions, it is easy to track these solutions numerically. We have carried out such an analysis and our results 
are displayed in Fig. 6 for r = 0.4084. The number of these coherent states increases as a function of K and r. We 
have also studied the stability of these states by carrying out a linear perturbation analysis of Eqs. (|3a[-^d|) around 
the phase locked solutions. The algebraic details of this analysis are given in Appendix. Numerical solution of the 
stability conditions shows that all these higher frequency states are stable (except for small bands of unstable regions 
which are indicated by dashed curves in Fig. 6). Thus these higher frequency states are genuine collective states of 
the system that are physically accessible. It should be mentioned that similar collective states were also observed 
in the delay coupled phase only model studied by Schuster et al. p2| ] . One difference between those results and our 
amplitude inclusive model is that the frequencies of our model are slightly lower due to amplitude effects. When the 
phase locked amplitudes of the two oscillators are equal or if their ratio is close to unity, the effect of the amplitude 
on the magnitude of the phase locked frequencies ceases. 



We end this analysis of the N — 2 model with a short description of the typical phase diagram for finite time 
delays. Such a diagram is shown in Fig. 7 for r = 0.4084 which reveals a much richer structure in comparison to 
the Aronson et al pl| diagram of Fig. 1. Note that one no longer has the clean separation of the Bar-Eli region, the 
phase locked region and the phase drift region into three disjoint regions that converge at a single degenerate point. 
Instead the phase locked region now always surrounds the Bar-Eli region and the single degenerate point is replaced 
by a series of X points resulting from the braided structure of the phase locked region in the vertical direction. The 
dotted curve (obtained numerically) marks the separation of the phase locked and the incoherent regions. This curve 
also represents the birth of two fixed points of the system, one stable and the other unstable. These branches can 
be seen in Fig. 8a where the phase locked amplitudes are plotted at A = 2. The dashed curves are the unstable 
branches and the solid curves stable branches. At large values of K other bifurcation curves appear in the phase 
locked region indicating the appearance of higher frequency states f3^] . Figs. 8b and 8c show respectively the phase 
locked solutions for A = 7 as K is varied and for K = 1 as A is increased. It should be noted that the basic nature of 
the transitions, namely a supercritical Hopf bifurcation, is preserved in the presence of time delay. Fig. 8b, a typical 
example, illustrates this clearly. 



III. ASSEMBLY OF N DELAY COUPLED OSCILLATORS 



A. Model Equations 



In this section we study the interaction of a large number of limit cycle oscillators that are globally coupled with a 
linear time delayed coupling. To describe such a system, we introduce the following set of model equations, 



G 



Z 3 (t) = (1 + iuj- | Z 3 (t) | a )^.(t) + f £ [3k (f - r) - Zj(t)] 

k=l 



N 



K' 
N 



(16) 



where j — 1, N, K' = IK is the coupling strength and r is the delay time. The last term on the R.H.S. has been 
introduced to remove the self-coupling term. The model is a straightforward generalization of the mean field model 
studied extensively by Ermentro ut p^j , Mirollo and Strogatz , and others and reduces exactly to their model 
for t = 0. Mirollo and Strogatz p3[ have provided rigorous analytical and numerical conditions for amplitude death 
in such a mean field model system. Their conclusions, in general, are similar to the case of TV = 2, namely, that one 
needs a sufficiently large variance in frequencies for death to occur and K has to be sufficiently large. In a short while 
we will discuss our findings for the time delayed model in the light of their results. 

As is customary in mean field models we can also define a centroid or 'order parameter' of the form, 



Z 



1 N 



(17) 



where R and <fi denote the amplitude and phase of the centroid. The order parameter is a useful quantity in the 
large N model since its behavior provides qualitative and quantitative clues about the collective and nonstationary 
states of the system, e.g. R = (in the large time limit) is indicative of an incoherent state whereas R = 1 marks a 
'phase locked' state. As has been noted by Matthews and Strogatz |2j||2(| the time behavior of R can also characterize 
chaotic states and other nonstationary states like large amplitude Hopf oscillations. We will also examine the behavior 
of this parameter in certain cases to track the effects of the time delay parameter on the collective dynamics of the 
large N system. In terms of the order parameter defined above the model equations can be written as follows. 



Zj(t) = (1 - K'd + itUj- | Zj(t) \ 2 )Z 3 (t) + K'Z(t - t) 



K[_ 
~~N 



Z 3 {t 



(18) 



where d = 1— 1/N. To study the stability of the origin, we once again carry out a linear perturbation analysis of (|l6| ) 
with the perturbation varying as e At . The linearized matrix of (jig) is given by 



B 



h I 
f h 

I f 



l-N 



(19) 



where l n = 1 — K'd + iuj n , / = jr e A ' 



That is 



B, 



1 - K'd 

N e ' 



if i 



3 



if* ^ j 



It is more convenient to cast the eigenvalue problem in terms of another matrix C where we define C = B + (K'd — 1)1 
(with / being the identity matrix). If fi is the eigenvalue of C then it is related to A by the relation, /i = \+{K'd— l).The 
matrix C is given by 



Cm. 7?. — 



det 



/ 



/ 



if m = n, 
if m 7^ n. 



The eigenvalues /i are obtained by solving det(C — [il) = 0, i 



ibJi- (jl f 
f iu 2 — [i 



(20) 



/ 



iujn — fi 



(21) 



This eigenvalue equation can be compactly expressed as a product of two factors, 



7 



JV 



JJ (iuj k — A* — /) 



1 + / 



Lfc=i 



AT 

E 



(22) 



As pointed out by Matthews and Strogatz |25|,g6|], the first factor represents the continuous spectrum of the system 
whereas the second factor gives us the discrete spectrum. In general it is not possible to solve the characteristic 
equation ( p2| ) analytically and for large N the numerical tracking of all the eigenvalues is also an arduous task. 
There are two interesting limits however, in which the analysis gets considerably simplified. If the N oscillators have 
identical frequencies then it is possible to obtain exact algebraic relations for the critical curves marking the amplitude 
death region. This is of considerable interest to us in view of the novel result of the N — 2 model which showed 
amplitude death for A = 0. The interesting question to ask is whether such a phenomenon exists for the arbitrary N 
case. We will study this question in the next subsection. Another interesting limit is the N — > oo limit, often called 
the thermodynamic limit. It is once again possible to obtain some exact analytic results in this limit rega rding the 
behavior of the critical curves in the bifurcation diagram. We will carry out this analysis in the subsection III C| . As 
shown by Mirollo and Strogatz |?3) (for the t = case) the infinite system results are of practical interest since they 
provide a fairly accurate picture of amplitude death for the large but finite system. A similar conclusion holds for the 
time delay case as well, as we find out by comparing the analytic results of the thermodynamic limit with numerical 
solutions for a large number of oscillators. In the final subsection we look at the time behavior of the order parameter 
in the various collective regimes and discuss its dependence on the time delay factor. We also briefly discuss the region 
of nonstationary solutions (chaos, Hopf oscillations etc. 
for the t — case. 



that was first pointed out by Matthews and Strogatz [E5_ 26 



B. Amplitude death of identical oscillators for arbitrary N 



In this section we treat the case of identical oscillators suffering amplitude death for an arbitrary number of globally 
coupled oscillators, N. For a set of identical oscillators the frequency distribution of the system is a delta function, 



g(uj) = 8(lo - loq) 



(23) 



where ujq is the natural frequency of each oscillator. With this assumption the matrix B becomes a circulant matrix 
whose eigenvalues can be expressed in terms of the N th roots of unity. Since in this particular case only two kinds of 
coefficients appear in the matrix B, the eigenvalues become much simpler. Explicitly they are given by 



\ = {1-K'd- 



UV 



K'de 



-At 



1 - K'd 



K Ati 

TV h 



(24) 



in which the second eigenvalue has a degeneracy of N — 1. Considering both the eigenvalue equations and following 
a similar procedure as described for the N = 2 case, we obtain the following set of critical curves: 



T a (n,K) 



n(n,K) 



2tvk + cos 1 [l — 

2(n + ljTr-cos- 1 [l- jfa] 



r c (n,K) = 



luq + ^2K'd 
2(n+ 1)tt 



1 



cos 



1-K'd 



K>(\- 



UJ 



T d (n,K) 



2mt 



d)} 2 - {K'd -If 



cos 



1-K'd 



K'{l-d) 



uj + ^/{K'(l-d)}Z-(K'd-iy 



(25a) 



(25b) 



(25c) 



(25d) 



Note that N enters as a parameter (through d = 1 — 1/N) in these sets of curves and the family of curves is now 
four instead of the two found for the N = 2 case. In fact for N = 2, the curves T a (n,K) and T c (n,K) combine to 
give Ti(n, K) and Tb(n, K) and Td(n, K) combine to give Ti(n, K) of the previous section. In Fig. 9 we display some 
typical death island regions for various values of N as obtained from the critical curves ( ^5a| )-(25d) with n = 0. The 
sizes of the islands are seen to vary as a function of N and approach a saturated size as N — > oo. The existence of 
these regions independently confirmed by direct numerical solution of the coupled oscillator equations, demonstrates 
that the amplitude death phenomenon for identical oscillators happens even in the case of an arbitrarily large number 
of oscillators. They also display multiple connectedness of the death region for higher values of u>o as was seen for the 
N = 2 case. 
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C. Oscillator death in the thermodynamic limit 



In the thermodynamic limit (N — ► oo, d — ► 0) the summation in the discrete spectrum of Eqn. (|22|) can be replaced 
by an integral and replacing /i by its definition, the discrete eigenvalue equation can be written as, 



— oo 



where g{u>) denotes the distribution of the frequencies. The continuous spectrum is given by 

A = 1 — K + iui where uj 6 9(&) (27) 

From (p7|) one can see that one of the critical curves is if' = 1. And the amplitude death region should lie to the 
right of this region in order for the eigenvalues to have negative real parts. To determine the critical curves from the 
discrete spectrum, we substitute a + i/3 for A in (|2^) , rationalize the denominator and equate the real and imaginary 
parts to zero, to obtain at a = 0, 

/oo j£l y 

■„ (K'-iy + ip-uy 9{uj) duj = l^ cos(/3t) ' (28a) 

/OO Li) Q 1 

These two relations define the other critical curve and for a given g(to) they provide a complete analytic description. 
As an example, let us choose a simple uniform distribution given by, 

g( U ) = /£' G . (29) 

v ' \0, otherwise v ; 

The integrals I\ and I2 can now be explicitly carried out and we get, 



tan 1 I " ,' I + tan 



(3 \ 1 / a - P \ y ( 2ff(K'-\ 



K'-lJ \K'-lJ \{K'-1) 2 

2 c 
K 



cosOSr), (30a) 



and 



Thus in the limit of an infinite number of oscillators the above set of transcendental equations, parametrized by the 
time delay, provide a complete description of the critical curve for amp litude death (along with K' = 1). It is also 
easy to see that for t = 0, j3 — is the solution one obtains from ( |30rj| ) and upon substitution of this in ( |30aD one 
obtains 

_i / cr \ <7 



tan "M^n - * ( 31 ) 



which is the equation for the critical curves in the no delay case as was obtained by |2j^,|23]]. In Fig. 10 we plot the 
region of amplitude death for a typical value of r. We also note that the general nature of the bifurcation curves is 
similar to the N = 2 case including the occurrence of amplitude death for A = 0. As mentioned earlier, Mirollo and 
Strogatz |23|] also point out that the infinite system provides a fairly accurate description of the large finite system. In 
order to check this out for the time delayed system we have carried out numerical studies of a large system (N = 800) 
of oscillators for a few values of K and A at a fixed value of r. The values of A and K at which the transition from 
the phase locked region to the death region occurs are marked with circles in Fig. 10. As we see these points lie 
nearly on the analytic curves for the N = 00 case and thus the bifurcation diagram for the finite large system can be 
expected to be similar to the infinite case. 
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D. Behavior of the order parameter 

It is also interesting to examine the behavior of the order parameter in the limit of large number of oscillators since 
it provides a good description of the macroscopic behavior of the system. Our particular interest now is to detect any 
signatures of the time delay effect that might show up in the temporal behavior of the order parameter. 

We begin by studying its characteristic behavior in the amplitude death region where it is known to spiral in time 
to the origin. To investigate the influence of the time delay we can measure this exponential decay rate for various 
values of r. Our numerical results are plotted in Fig. 11 where the decay exponent is plotted against r for two values 
of K and for a particular value of co for the case of N = 800 oscillators. The range of r spans the width of a death 
island. The exponent is seen to have a maximum negative value in the centre of the island and it gets less negative 
as one moves towards the boundaries. Thus time delay has a dynamical influence on the rate at which the oscillators 
collapse into the origin with the maximum rate of decay occurring in the centre of the island. 

The death islands are surrounded by the phase locked regions in which the order parameter assumes a non zero 
constant value. The magnitude of this value is found to be a function of r. Our numerical results show that the 
saturated value of the order parameter grows in an algebraic fashion as a function of r as we move from the death 
island boundary into the phase locked region. This is demonstrated in Fig. 12 for K' = 5. 

Finally we examine the behavior of the order parameter in the narrow region straddling the boundary between the 
phase locked and the phase drift region where Matthews and Strogatz |^5|,^6| have found non-stationary behavior 
corresponding to Hopf oscillations, large oscillations, quasiperiodicity and chaos. We find that the order parameter 
displays all these characteristic behaviors even in the presence of the time delay. However there is an overall reduction 
of the phase space area of this region compared to the no delay case. To see this clearly we choose to consider the 
temporal record of the amplitude of the order parameter. The value of the order parameter is constant in the phase 
locked region, and tends to zero as l/y/~N with N, the number of the oscillators, in the incoherent region. In the 
irregular region a variety of behaviours can be seen. We take a long time record of R discarding the first few hundred 
points. The average value of this record, < R >, as shown in Fig. 13, for N = 500, shows a plateau in the irregular 
region. This is almost unchanged when the number of oscillators is increased from N — 500 to N = 5000. This region 
shrinks considerably as the delay parameter r is increased. We have also observed other qualitative changes in the 
nature of the time variation of the order parameter which suggest subtle influences of time delay on the dynamical 
pattern of the nonstationary states. A more detailed and quantitative study of these effects are in progress and will 
be reported elsewhere. 



IV. SUMMARY AND DISCUSSION 

We have studied the effect of time delay on the collective dynamics of a system of coupled limit cycle oscillators that 
are close to a supercritical Hopf bifurcation. The oscillators are linearly coupled in a global fashion and time delay is 
introduced in the coupling term. The principal effects of time delay on the stationary collective states like amplitude 
death and phase locking are clearly evidenced even in a simple two coupled oscillator system. The phase diagram 
shown in Fig. 7 summarizes these results. The boundary curve for the amplitude death region is a sensitive function 
of the time delay parameter and the phase space distribution of locked states, amplitude death and incoherent states 
is quite complex. The phase locked region shows the existence of multiple locked states all of which are stable. A 
surprising result is the occurrence of amplitude death even for two identical oscillators over a range of time delay 
values (Fig. 5). Such death islands in the parametric space of the coupling strength and the time delay factor are 
found to persist even for a large number of identical oscillators. Our numerical findings are backed by analytic results 
of linear stability of the amplitude death and phase locked states for the N — 2 and the N — ► oo cases. The large N 
results are found to agree very closely with the N — > oo curves. For the case of identical oscillators we provide exact 
analytic curves for the island boundaries. Our principal focus has been on the investigation of the stationary states 
but we have also confirmed the existence of the nonstationary states (e.g. Hopf oscillations, chaos) in the finite time 
delay model. We hope to carry out a more detailed investigation of these states in the future. 

As we have mentioned earlier, time delay effects in the context of the coupled oscillator model have not been looked 
at before. In view of the vast number of applications of the coupled oscillator model our results may have significance 
for some biological or physical systems. There are many physical examples of amplitude death in real systems. One of 
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the earliest that was investigated both theoretically and experimentally is that of coupled chemical oscillator systems 
e.g. coupled Belousov-Zhabotinskii reactions carried out in coupled stirred tank reactors They can also occur 

in ecological contexts where one can imagine two sites each having the same predator-prey mechanism which causes 
the number density of the species to oscillate. If the species are capable of moving from site to site at a proper rate 
(appropriate coupling strength) the two sites may become stable (stop oscillating) and acquire constant populations. 
Another important application of this concept is in pathologies of biological oscillator networks e.g. an assembly of 



cardiac pacemaker cells 14 . Amplitude death signifies cessation of rhythmicity in such a system which is otherwise 
normally spontaneously rhythmic for other choices of parameters. For the onset of such an arrhythmia, current models 
based on coupled oscillator networks need to assume a significant spread in the natural frequencies of the constituent 
cells (oscillators) ^3|. Our result of amplitude death for identical oscillators demonstrates that this assumption may 
not be necessary if one takes into account time delay effects arising naturally from the finite propagation times of 
the signals exchanged between the cells. Another possible application is in the area of high power microwave sources 
where it is proposed to enhance the microwave power production by phase locking a large number of sources such 
as relativistic magnetrons ||. Time delay effects, arising from the finite propagation time of information signals 
traveling through the connecting waveguide bridges, could impose important limitations on the connector lengths and 
geometries in these schemes. Our findings could provide a guideline in this direction. Coupled oscillators also feature 
in several neural network configurations where multiple phase locked states arising from time delay effects could play 
an important role. 



Finally we would like to discuss some further extensions of the present work and possible future directions of 
research in this area. Our present results have been largely derived from the simple model of limit cycle oscillators, 
which are close to a supercritical Hopf bifurcation and which are linearly coupled in a global fashion. It would be 
interesting to carry out a similar analysis for limit cycle oscillators with local coupling. For the particular case of death 
of identical oscillators we have confirmed that the result holds for the locally coupled model as well but a more general 
investigation of other collective states remains to be done. Likewise the introduction of more complicated dynamics 
in the individual oscillators, such as shear (amplitude dependent frequency) [fl9f and higher order nonlinearities can 
open up rich possibilities for the interplay of time delay and the collective dynamics of the system. 



APPENDIX A: STABILITY OF THE FREQUENCY LOCKED SOLUTIONS 



The phase locked solutions for the case of N = 2 are given by (ipi, tp 2 , if>3, ipi) = (Ri, R 2 ,Qt + cti,flt + a 2 ). To 
determine the stability of the phase locked states we carry out a linear perturbation analysis about these solutions. 
Suppose A represents the eigenvalue. The characteristic equation is then given by 
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where the coefficient functions are given by 

/oo = ABC X C 2 + D\D\ - AC 2 DlR x /R 2 - BC 1 D 2 1 R 2 /R 1 , 
/01 = -Add - BC\C 2 - AD\R\jR\ - BD\R%/R\ 

+ ABC 1 R l /R 2 + ABC 2 R 2 /R l + C 2 D 2 R 1 /R 2 + dDfR 2 /Ri, 
/02 = AB + dC 2 + D 2 R 2 /R 2 + D 2 R 2 /R 2 - AC l R 1 jR 2 - BC 1 R 1 /R 2 
— AC 2 R 2 j R\ — BC 2 R 2 j Ri, 

A-B + C 1 R 1 /R 2 + C2R2/R1, 







/o4 
/20 



-ABC X C 2 - C 2 C 2 - dC 2 D 2 - ddD 2 
C\D Y D 2 - ClD x D 2 + 2ddDiD 2 - 2D\D\ 

AdD 1 D 2 R 1 /R 2 + AC 2 D 1 D 2 Ri/R2 + BdD 1 D 2 R 2 jR 1 - BC 2 D 1 D 2 R 2 /R 1 



(Al) 



where A = 1-K-3R(, B = 1-K — 3R 2 , C x = K cos(-Or + a), C 2 = Kcos(Q,T+a), D l = K s\n(-VLr + a) , D 2 = 
K cos(— fir — a). We can also write the above determinant in the form of the following equation: 
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AC x D\R x jR 2 + BC 2 D\R 2 /R X , 
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C\C 2 R\j R 2 — C\C\R 2 f R\ 



+ C 1 D 1 D 2 R 1 /R 2 + C 2 D X D 2 R 2 IR X 
- 2C 2 D 2 1 R 2 /R 1 - 2C 1 D 1 D 2 R 2 /R 1 , 
= -2C 1 C 2 + 2D 1 D 2 , 



2C 2 D 1 D 2 R 1 /R 2 - 2C 1 DjR 1 /R 2 
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The above equation (A2) has been solved numerically to generate the stability curves displayed in the paper. 
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FIG. 1. Bifucrcation diagram of Aronson, et al. of two coupled oscillators without time delay in the plane of coupling 
strength (K) and the frequency mismatch (A). Region I, bounded by the curves K = 1 and K = (1 + A /4)/2, represents the 
oscillator death region. Region II, bounded by K — A/2 (when A < 2) and K = (1 + A 2 /4)/2 (when A > 2), represents the 
frequency locked state. Region III represents the incoherent or drift solutions. All three regions meet at a degenerate point 
(Jf,A) = (l,2). 
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FIG. 2. The effect of time delay on the boundary of the amplitude death region of two coupled oscillators. The critical 
curves are plotted from Eqs. (paHsd) for Hi = 10. Region I represents the amplitude death region. For small values of r the 
denerate point (K, A) = (1, 2) disappears and the death region is bounded by the curves S- (solid lines) and S+ (dashed lines) 
as defined in the text. As the delay parameter r is increased, the bounding curves get deformed continuously. 
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FIG. 4. The existence of multiple death islands. In the (r, K) space the amplitude death region is multiply connected (i.e. 
higher order death islands exist) for higher values of ui . The figure shows that for ui — 30, there are three death islands which 
are defined by the set of curves Ti^im, K) where m = 0, 1, 2. 
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FIG. 5. Dependence of the size of the death island on the common intrinsic frequency, loq. Below a certain threshold, which 
is given by the condition of intersection of the curves n(0, K) and T2(0, K), the amplitude death region disappears. This value, 
found numerically, is 4.812. 
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FIG. 6. The phase locked solutions. The common frequency Q, the phase difference a and the amplitudes of the individual 
oscillators plotted as a function of the coupling strength, K for r = 0.4084, A = 1, and u> = 10. The oscillators have the phase 
difference a, either around (i.e. in-phase synchronization) or around 7r (anti-synchronization). The unstable portions of the 
solutions are plotted as dashed lines. 




FIG. 7. The bifurcation diagram in the (K, A) space for r = 0.4084 and Q = 10. The region marked I is the amplitude death 
region, II corresponds to the phase locked solutions, and the region III is the phase drift or incoherent region. The dotted 
curve which separates the phase locked region from the incoherent region is obtained by numerical integration of the original 
equations. Notice that the degenerate points (marked by stars) have reappeared for this value of the delay parameter r. At 
higher values of K, other bifurcation curves appear indicating higher order frequency states as shown in Fig. 6. 
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FIG. 8. The amplitude curves _Ri and R2 of the phase locked solutions of the two oscillators plotted in three different regions. 
Ri corresponds to the amplitude of the oscillator with the smaller of the intrinsic frequencies and R2 that of the larger of the 
two. (a) At A = 2, the stable branches (solid curves) continue to exist whereas the unstable (dashed curves) of the two 
oscillators merge with the origin as K is increased, (b) The amplitudes of the oscillators in one of the 'closed loops' in the 
bifurcation diagram (previous figure) at A = 7. The origin is unstbale (dashed line) in the region where the periodic orbits are 
stable and stable outside this region, illustrating the supercritical Hopf bifurcation. And (c) shows the phase locked solutions 
as one moves along the line K = 1. In the first gap there are phase drift solutions and in the second gap amplitude death 
solutions. 
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FIG. 9. The existence of death islands for an arbitrary nu mber of globally coupled oscillators. The figure shows the death 
island regions for N — 2, 3, 4 and 00 plotted from the Eqs. (25a -25d) with n — 0. All the oscillators are assumed to have an 
intrinsic frequency of ljq = 10. The death island survives even in the limit of N — > 00. 
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FIG. 10. The amplitude death region for N globally coupled oscillators, as N 
frequencies, g(uj), with the average frequency, G) = 10. plotted from the Eq. (30a) and (30b) by eliminating (3 for r 
The death region is bounded by K' = 1 and taxi{a / K 1 ) = <r/(K' — 1) f° r t = 0. Both the curves meet at a degenerate point 
(K',cr) — (l,7r/2). As the delay parameter r is increased the boundary K' — 1 (i.e. K = 0.5) still continues to be one of the 
boundaries. But the second curve distorts in a continuous fashion and at a certain value of r it touches the a = axis, showing 
the death of identical oscillators. The points marked indicate the numerically found value of the boundary for N = 800. 



N = 800, a = 0,© = 10 




FIG. 11. The rate of approach to the origin of the amplitude of the order parameter, R. If we assume that R oc e 



V(r)t 



the 



exponent T]{t) plotted in the above graph shows that the rate of approach to the origin is maximum in the middle of the death 
island. As one approaches the boundaries of the island on either side the rate decreases to 0. 
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FIG. 12. Power law behavior of the magnitude of the order parameter at the boundary of the death island. The order 
parameter is inside the death island region due to the fact that the amplitudes of all the oscillators are zero. The order 
parameter shows a power law behavior at the transition between the death island and the phase locked region. The curve is 
approximated numerically, at K' = 5, by r = 0.432 + 0.0992(i? 3 + 0.022.R 2 + 0.154R). 
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FIG. 13. The effect of the time delay, r, on the irregular region which was found for r = by Matthews and Strogatz (see 
text). This region which appears as a plateau for r = as a whole shrinks as r is increased. The order parameter is unstable 
in this region. However a long temporal record of R gives nearly constant average value. The average value distinguishes the 
three distinct regions quite accurately. To the right of this plateau region is the phase locked region, and to the left is the 
incoherent region. 
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